Problema 1¶
Considere una partícula constreñida a moverse en un círculo de radio R. De esta manera la energía del sistema es $E=\frac{p^2}{2mR^2}$, lo cual corresponde a un movimiento con velocidad angular constante $\omega=\dot{\theta}=\frac{p}{mR^2}$. Ahora, consideremos que el sistema es perturbado periódicamente cada tiempo $\tau$ por una fuerza impulsiva. Esta fuerza tiene la misma intensidad $\epsilon$ y siempre apunta en la misma dirección y sentido. Esto significa que la fuerza en dirección tangencial dependerá de la posición de la partícula en el tiempo en el que se aplique.
Esta perturbación tiene el efecto de cambiar el momento angular de manera instantánea. Si la partícula se encuentra en la posición $\theta(t)$, entonces habrá un cambio en la cantidad de movimiento $\Delta p=\epsilon cos\theta$. Entonces, al iniciar el movimento, el momento permanecerá constante un tiempo $\tau$ y entonces tendrá un cambio instantáneo $\Delta p$ y permanerecá constante hasta que se cumpla otro periodo. Además, dentro del intervalo de tiempo donde no hay perturbación, la partícula se mueve conforme a $\theta(t)=\theta_0+\omega t$.
Para anializar este sistema se utiliza un método similar al mapeo de Poincare. Dado que la perturbación tiene un periodo $\tau$, podemos tomar un punto en el espacio fase y mapearlo al punto en el que estará un tiempo $\tau$ después.
Considere que las coordenadas iniciales antes de la perturbación son $\theta_0,p_0$ y escoja un valor de $\tau$.
a) Muestre que en el mapeo, los valores subsecuentes son:
$$\theta_{i+1}=\theta_i + \frac{p_{i+1}}{m R^2} \tau \quad \text{(módulo}\ 2\pi), \quad y \quad p_{i+1}=p_i + \epsilon cos\theta_i,$$donde se ha tomado en cuenta que la fuerza actúa cuando la posición es $\theta=\theta_i$. En adelante considere que $m$=1 y $R$=1.
b) Con $\epsilon=0$ utilize 200 condiciones iniciales arbitrarieas tales que $p_0 \in[-\pi/2,\pi/2]$ y $\theta_0 \in [0,2\pi]$ y encuentre las curvas en el plano $(\theta,p)$ que resultan de hacer muchas iteraciones del mapeo. Explique el resultado del mapeo y grafíquelo, tome por lo menos 2,000 iteraciones para cada condición inicial.
c) ¿El sistema tiene puntos fijos para $\epsilon \ne 0$?. De ser así, encuentrelos y explique el significado de dichos puntos. Los puntos fijos son aquellos que cumplen $\theta_{i+1}=\theta_i$ y $p_{i+1}=p_i$.
d) Realice simulaciones con $\epsilon = 0.01,0.05,0.1,0.2,0.3,0.4,0.7,1.1$. Explique detalladamente sus resultados.
Solución¶
Observamos que la energía en este caso es exclusivamente cinética. Pero, además, vemos que la energía dada por el problema tiene un factor de $\frac{1}{R^2}$ adicional. Por claridad, haremos el desarrollo para obtener tanto las ecuaciones de movimiento como los momentos en coordenadas polares $(r,\theta)$.
En primer lugar, sabemos que el potencial en este caso es $V(r,\theta)\equiv 0$. Con esto, tanto el lagrangiano $\mathcal{L}$ como el hamiltoniano $\mathcal{H}$ (o la energía) del sistema se reducen símplemente al término de la energía cinética. La única diferencia será que el primero está expresado en términos de las posiciones y las velocidades, mientras que el segundo está escrito en términos de las posiciones y momentos generalizados. Es decir:
$$\begin{cases}\mathcal{L}(r,\theta,\dot{r},\dot{\theta}) &=& T-V &=& T &=& \frac{mv^2}{2} \\ \mathcal{H}(r,\theta,p_r,p_{\theta})&=& T+V &=& T &=&\frac{p^2}{2m} \end{cases}$$.
Ahora, en la tarea examen anterior habíamos visto el desarrollo para la deducción de la velocidad y aceleración en coordenadas polares. Vimos que, si tomamos que $\vec{r}(t)=r(t)\hat{r}(t)$,con $\hat{r}=(cos\theta(t),sen\theta(t))$, entonces la velocidad y la aceleración son: $\begin{cases}\vec{v}=\dot{\vec{r}}=\dot{r}\ \hat{r} + r\dot{\theta}\ \hat{\theta} \\ \vec{a}=\ddot{\vec{r}}=(\ddot{r}-r\dot{\theta}^2)\ \hat{r}+(2\dot{r}\dot{\theta}+r\ddot{\theta})\ \hat{\theta} \end{cases}$, en donde tomamos $\hat{\theta}=(-sen\theta(t),cos\theta(t))$.
De la expresión de la velocidad vemos que, en coordenadas polares se tiene que $v^2=\dot{r}^2+r^2\dot{\theta}^2$, por lo que $\mathcal{L}=T=\frac{m}{2}(\dot{r}^2+r^2\dot{\theta}^2)$. Entonces, vemos que:
$$\begin{cases}p_r &=& \frac{\partial \mathcal{L}}{\partial \dot{r}} &=& m\dot{r} \\ p_{\theta} &=& \frac{\partial \mathcal{L}}{\partial \dot{\theta}} &=& m r^2\dot{\theta} \end{cases} \quad .$$Entonces, al despejar $\dot{r}$ y $\dot{\theta}$ de estas expresiones y sustituir en el lagrangiano se obtiene que $\mathcal{H}=T=\frac{p_r^2}{2m}+\frac{p_{\theta}^2}{2mr^2}$. Entonces, $$E=\mathcal{H}=\frac{r^2 p_r^2 + p_{\theta}^2}{2mr^2}.$$
Esto explica el factor de $1/R^2$ que aparecía originalmente, sin embargo, preferiremos trabajar con esta expresión para la energía y su análogo en el lagrangiano.
A continuación vemos que, la constricción que tenemos está dada por $f(r,\theta)=r-R=0$ . Dado que las ecuaciones de movimiento las podemos obtener mediante $\frac{d}{dt}\big(\frac{\partial \mathcal{L}}{\partial \dot{q}_i} \big)=\frac{\partial\mathcal{L}}{\partial q_i}+ \lambda \frac{\partial f}{\partial q_i}$, obtenemos que las ecuaciones de movimiento son: $\begin{cases} \frac{d}{dt}(mr^2\dot{\theta}) &=& 0\\ \frac{d}{dt}(m\dot{r}) &=& mr\dot{\theta}^2 + \lambda \end{cases}$. Al utilizar la constricción se obtiene, como era de esperarse que:
$$\begin{cases}r &=& R \\ \theta &=& \theta_0 + \omega t \\ \lambda &=& -mR\omega^2 \end{cases} \quad .$$Esto nos dice que en el espacio fase, mientras no se considere la perturbación, se tiene que: $\begin{cases} r=R=cte \\ p_r = m\dot{r} = 0 \\ \theta=\theta_0 + \omega t \\ p_{\theta} =mr^2\dot{\theta}=mR^2\omega = cte \end{cases}$. Entonces, sólo tendrá relevancia fijarse en el plano $(\theta,p_{\theta})$ para este problema.
a) Ahora consideraremos los fuerza impulsiva en el problema. En principio esta es de la forma $\vec{f}=\kappa \hat{\mu}$, donde $\hat{\mu}=(cos\alpha,sin\alpha)$ es un vector unitario en dirección arbitraria dada por $\alpha$. Ahora, debido a que la fuerza de constricción $\lambda \frac{\partial f_{ctr}}{\partial r}=-mR\omega^2$ se ejerce en dirección radial, la componente radial de la fuerza impulsiva $\vec{f}$ no afectará el movimiento de la partícula. Entonces, la componente en dirección $\hat{\theta}$ de la fuerza esta dada por $\vec{f} \cdot \hat{\theta}=\kappa(-cos\alpha sin \theta + sin\alpha cos \theta)$. En particular, si $\alpha=\pi/2$, entonces $f_\theta=\vec{f}\cdot \hat{\theta}=\kappa cos\theta$. Sabemos que $f_\theta=\frac{\Delta p_\theta}{\Delta t}$, i.e. $\Delta p = f \Delta t$. Como el cambio del momento $\Delta p_\theta$ es finito y $\Delta t$ es infinitesimalmente chico, entonces la fuerza deberá ser infinita para darnos un cambio en el momento de la forma $\Delta p_\theta = \epsilon(-cos\alpha sin \theta + sin\alpha cos \theta)$. Si $\alpha=\pi/2$ se tiene el caso mencionado en el problema: $\Delta p_\theta =\epsilon cos \theta$.
Ahora, consideremos que la partícula tiene posición $(\theta_i,p_{\theta i})$ en el espacio fase en el momento en el que se le aplica la fuerza. Entonces, su momento después de la interacción es $p_{\theta i+1}=p_{\theta i} + \Delta p_\theta=p_{\theta i} + \epsilon(cos\alpha sin \theta + sin\alpha cos \theta)$. Además, como $p_{\theta}=mR^2\omega$ y $\theta=\theta_0+\omega t$, al despejar $\omega$ y sustituirla en la expresión de $\theta$ se tiene que $\theta= \theta_0 + \frac{p_\theta}{mR^2}\tau$.
Entonces, el mapeo del punto $\theta_i$ al siguiente punto en el que se aplica la fuerza $\theta_{i+1}$ será:
$$\theta_{i+1}=\theta_i + \frac{p_{\theta i+1}}{m R^2} \tau,\quad \text{donde} \quad p_{\theta i+1}=p_{\theta i} + \epsilon (-cos\alpha sin \theta + sin\alpha cos \theta).$$Como ya mencionamos, si $\alpha=\pi/2$, entonces $p_{\theta i+1}=p_{\theta i} + \epsilon cos \theta$.
b) A partir de aquí consideraremos m=1,R=1. Ahora, dado valores arbitrarios de $\tau$ y $\alpha$ haremos 2,000 iteraciones a partir de una condición inicial $(\theta_0,p_{\theta 0})$ en el espacio fase. Esto se hará para 200 condiciones iniciales tales que $p_0 \in[-\pi/2,\pi/2]$ y $\theta_0 \in [0,2\pi]$.
Para empezar, definimos una función que, dada la condicón inicial nos regrese un array con todos los valores de las distintas iteraciones del mapeo.
import numpy as np
import pylab as pl
import math as mt
#from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
%matplotlib inline
def mapeo_iterado(theta_0,p_0,iteraciones=2000,m=1.,R=1.,tau=5.,epsilon=0.,alpha=np.pi/2):
theta=np.zeros(iteraciones)
p_theta=np.zeros(iteraciones)
if theta_0<0:
theta_0+=2*np.pi
theta[0]=mt.fmod(theta_0,2*np.pi)
p_theta[0]=p_0
for i in xrange(iteraciones-1):
p_theta[i+1]=p_theta[i]+epsilon*(-np.cos(alpha)*np.sin(theta[i])+np.sin(alpha)*np.cos(theta[i]))
theta[i+1]=theta[i]+p_theta[i+1]*tau/(m*R**2)
theta[i+1]=mt.fmod(theta[i+1],2*np.pi) #tomamos theta modulo 2pi (puede ser negativo)
if theta[i+1]<0:
theta[i+1]+=2*np.pi #volvemos positivo el ángulo
return theta,p_theta
Ahora tomamos una serie de condiciones iniciales. Variamos 20 valores del ángulo y 11 valores de momento para cada uno de ellos, por lo que en total son 220 condiciones iniciales, del lado derecho de plano pues se tomaron sólo ángulos positivos.
angulos=np.linspace(0,2*np.pi,20)
p_angulos=np.linspace(-np.pi/2,np.pi/2,11)
sols_theta=[]
sols_pth=[]
fig=plt.figure(figsize=(12,12))
for i in xrange(len(angulos)):
for j in xrange(len(p_angulos)):
theta0=angulos[i]
pth_0=p_angulos[j]
th,p=mapeo_iterado(theta0,pth_0,iteraciones=2000)
sols_theta.append(th)
sols_pth.append(p)
plt.plot(th[0],p[0],marker='o')
plt.plot(th[-1],p[-1],marker='x')
plt.plot(th,p)
if min(p)!=max(p):
print "Cambio de Momento" #verificamos que no hay cambio en el momento puesto que epsilon=0
th0,p0=mapeo_iterado(angulos[0],p_angulos[0],iteraciones=2000)
plt.plot(th0[0],p0[0],marker='o',color="black",label="Pto. Inicial")
plt.plot(th0[-1],p0[-1],marker='x',color="black",label="Pto. Final")
plt.ylabel('$p_{\\theta}$',fontsize=20)
plt.xlabel('$\\theta$',fontsize=20)
plt.legend()
plt.grid()
plt.show()
En este caso observamos que nunca hay cambios en el momento debido a que $\epsilon=0$. Esto nos dice que el movimiento es circular uniforme durante todo el tiempo. Por tanto, las curvas en el plano $(\theta,p_\theta)$ que resultan de aplicar iteradamente el mapeo son líneas horizontales entre $0$ y $2\pi$. El hecho de que un mapeo iterado pueda cubrir por completo una linea horizontal o no depende del número de iteraciones, el valor del intervalo $\tau$ y del valor del momento inicial. Por ejemplo, observamos que hay 3 rectas donde no se dibujaron líneas. Estas son para $p_{\theta\ 0}=-\frac{2\pi}{5},0,\frac{2\pi}{5}$. En el caso de $p_{\theta\ 0}=0$ es evidente que no hay cambio debido a que el momento inicial es cero. Sin embargo, para los otros casos observamos que $\tau p_{\theta\ 0}=\pm 2\pi$, por lo que en cada iteración la partícula da exactamente una vuelta y regresa a su posición inicial.
c)El mapeo puede separarse en tres términos: $\theta_{i+1}=\theta_i+\frac{p_{\theta i}}{mR^2}\tau + \frac{\epsilon}{mR^2}\tau(-cos\alpha sin \theta + sin\alpha cos \theta)$, donde los dos primeros corresponden a movimiento circular uniforme. El tercer término corresponde a la perturbación periódica, donde consideramos un valor arbitrario de $\epsilon$ distinto de cero.
Para que $p_{\theta i}$ este fijo en el mapeo debe pasar que el tercer término sea cero (en vez congruente con $2\pi$), i.e. debe pasar que $\theta_i=\pm \alpha$. Esto significa que la fuerza impulsiva tiene solo dirección radial y por tanto no contribuye en un cambio de momento en dirección tangencial.
Con esto, para que $\theta_i$ sea un punto fijo en el mapeo, basta que sea un punto fijo del movimiento sin perturbación, que analizamos en el inciso anterior.
Por tanto, los puntos fijos serán aquellos en los que $\theta_i=\pm \alpha$ y $p_{\theta i}\tau \equiv 0\ (mod 2\pi)$. Notemos que esto no depende del valor de $\epsilon$.
d) Ahora tomaremos una serie de mapeos distintos con $\epsilon = 0.01,0.05,0.1,0.2,0.3,0.4,0.7,1.1$. Haremos una gráfica para cada uno de ellos variando las condiciones iniciales de modo que $p_0 \in[-\pi/2,\pi/2]$ y $\theta_0 \in [0,2\pi]$. Sin embargo, para evitar saturar la gráfica tomaremos solamente 45 condiciones iniciales y 20 iteraciones.
Para fines ilustrativos, deseamos graficar las trayectorias así como sus correspondientes puntos finales e iniciales. Para que estos tengan el mismo color para una de condición inicial dada, utilizaremos la función set_color_cycle, iterando sobre un array que repita 3 veces cada valor del código de colores.
epsilons=[0.01,0.05,0.1,0.2,0.3,0.4,0.7,1.1]
angulos=np.linspace(0,2*np.pi,9)
p_angulos=np.linspace(-np.pi/2,np.pi/2,5)
fig=plt.figure(figsize=(20,35))
fig.suptitle('Espacio Fase con Distintas $\\epsilon$',fontsize=18)
#plt.subplots_adjust(wspace=0.2)
#plt.subplots_adjust(hspace=0.3)
for i,eps in enumerate(epsilons):
ax = fig.add_subplot(4,2,i+1) #haremos un arreglo de 4x2 subplots
colormap = plt.cm.prism
colors=np.linspace(0,1.,45) #rango de colores para 5*5=25 condiciones iniciales
iteracion=np.zeros(3*len(colors))
for k in xrange(len(colors)):
iteracion[3*k:3*k+3]=np.array([colors[k],colors[k],colors[k]]) #triplicamos cada entrada del rango de colores
plt.gca().set_color_cycle([colormap(j) for j in iteracion]) #definimos la iteracion del ciclo de colores
for angini in angulos:
for pini in p_angulos:
th0=angini
p0=pini
th,p=mapeo_iterado(th0,p0,iteraciones=20,m=1.,R=1.,tau=5.,epsilon=eps,alpha=np.pi/2)
plt.plot(th,p)
plt.plot(th[0],p[0],marker='o')
plt.plot(th[-1],p[-1],marker='x')
plt.grid()
plt.ylabel('$p_{\\theta}$',fontsize=20)
plt.xlabel('$\\theta$',fontsize=20)
ax.set_title('$\\epsilon$ ={}'.format(eps),fontsize = 16)
plt.show()
Recordamos que el mapeo puede escribirse como: $\theta_{i+1}=\theta_i+\frac{p_{\theta i}}{mR^2}\tau + \frac{\epsilon}{mR^2}\tau(-cos\alpha sin \theta + sin\alpha cos \theta)$, por lo que la dinámica del sistema dada por la fuerza impulsiva depende del producto $\epsilon \tau$ que controlamos variando únicamente el valor de $\epsilon$.
Observamos que, para valores pequeños de $\epsilon$, el sistema se comporta de forma similar que en el inciso b), en donde las órbitas se mantienen horizontales. Esto significa que la fuerza tiene poco efecto en cambiar su momento angular. Conforme se aumeta la intensidad de la fuerza, el cambio en el momento angular es más significativo y las trayectorias comienzan a abarcar un mayor rango vertical en el espacio.
En este caso, tomamos $\alpha=\frac{\pi}{2}$, lo cual nos dice que la fuerza impulsiva es vertical hacia arriba. Además, dado que esta se aplica en la posición actual de la partícula, esta puede tanto aumentar como disminuir el momento angular en la trayectoria. Esto hace que las trayectorias puedan subir y bajar en el espacio conforme se avanza en las iteraciones, lo cual se puede observar en los mapeos con valores grandes de $\epsilon$. En estos casos, también ocurre que los valores del momento angular se dispersan fuera de la región de las condiciones iniciales en $p_0 \in[-\pi/2,\pi/2]$.
Observamos que hay un punto fijo en la condición inical $(\theta,p_\theta)=(\frac{\pi}{2},0)$ y $(\theta,p_\theta)=(\frac{3\pi}{2},0)$ , lo cual se debe a que la fuerza tiene dirección radia, como se discutió en el inciso c).
Se puede ver que se forma un pequeño círculo para la condición inicial $p_\theta=0$ y valores bajos de $\epsilon$ (en puntos que no son fijos) que se rompe en las útlimas dos gráficas. Esta pareceser una órbita cerrada en el mapeo debida su periodicidad. El hecho que se rompa para mayores valores de $\epsilon$ nos dice que el sistema se vuelve más caótico conforme aumenta este parámetro y se debe a que se rompe la periodicidad en el movimiento de la partícula alrededor del círculo.
Problema 2¶
Considere un sistema dado por $H=\frac{1}{2}(p_x^2+p_y^2)+\frac{1}{2}(x^2+y^2)+x^2y-\frac{1}{3}y^3$.
En este caso, la energía $E=H$ es una constante de movimiento.
a)Muestre que el sistema tiene órbitas acotadas. ¿Para qué energías se dan dichas órbitas? De aquí en adelante consideremos condiciones iniciales con energía $E<\frac{1}{6}$.
b) Realice un programa que elija condiciones iniciales con $x=0$ y valores de $(y,p_y)$ compatibles con la condición de energía dada en el inciso a). Escoja un valor para la energía y ella determinará un valor de $p_x$ para esta condición inicial. Ahora obtenga la trayectoria para el valor inicial que se escoja y resuelva las ecuaciones de movimiento dadas por: $$\dot{x}=\frac{\partial H}{\partial p_x},\quad \dot{y}=\frac{\partial H}{\partial p_y},\quad \dot{p_x}=-\frac{\partial H}{\partial x}\quad \text{y} \quad \dot{p_y}=-\frac{\partial H}{\partial y}$$
c)Ahora programe una función que obtenga todos los puntos de la trayectoria solución cuando $x=0$. Escoja 15 trayectorias tales que $x=0$, $y \in[0,0.1]$ y $p_x \in[-0.15,0.15]$. El valor de $p_y$ deberá se compatible con la condición de energía del siguiente inciso.
d) Realice simulaciones con las condiciones iniciales anteriores para los siguientes valores de energía $E=0.01,0.03,0.1,0.12,0.15,0.166$. Haga la gráfica del espacio fase $(y,p_y)$ para cada condicón inicial. De esta manera, podemos obtener el comportamiento del sistema cuando barremos valores de la energía que se van acercando a $1/6$. Explique los resultados obtenidos.
Solución¶
a)Para resolver este problema integraremos las ecuaciones de movimimiento en un intervalo de tiempo largo y veremos si estas estan contenidas dentro de una circunferencia de radio R en el espacio de posiciones, viendo en cada caso la energía que estas tienen.
Las ecuaciones de movimiento estan dadas por: $\quad \begin{cases} \dot{x}=\frac{\partial H}{\partial p_x} &=& p_x \\ \dot{y}=\frac{\partial H}{\partial p_y} &=& p_y\\ \dot{p_x}=-\frac{\partial H}{\partial x}&=& -(x\ +\ 2xy) \\ \dot{p_y}=-\frac{\partial H}{\partial y}&=& -(x^2\ + \ y\ - \ y^2) \end{cases}$. A continuación programaremos las ecuaciones de movimiento, la función de energía y la norma de los puntos de la trayectoria:
from random import *
import scipy.integrate as spi
def ec_mov(x,t): #x=(x1,x2,x3,x4)=(x,y,px,py) asumiremos que este es un array
dx=x[2]
dy=x[3]
dpx=-(x[0]+2*x[0]*x[1])
dpy=-(x[0]**2+x[1]-x[1]**2)
return np.array([dx,dy,dpx,dpy]) #regresa un array
def funcion_energia(pini): #la entrada sera un array
x=pini[0]
y=pini[1]
px=pini[2]
py=pini[3]
energia=0.5*(px**2+py**2)+0.5*(x**2+y**2)+y*x**2-(1/3.)*y**3
return energia
def norma(x):
return np.sqrt(x[:,0]**2+x[:,1]**2)
Ahora, graficamos el potencial (los términos que depende de x,y):
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
fig = plt.figure(figsize=(15,15))
ax = fig.gca(projection='3d')
X = np.arange(-1, 1, 0.05)
Y = np.arange(-1, 1, 0.05)
X, Y = np.meshgrid(X, Y)
#Z = 0.5*(X**2+Y**2)+Y*X**2-(1/3.)*Y**3
Z =0.5*(X**2+Y**2)+Y*X**2-(1/3.)*Y**3
surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.coolwarm,linewidth=0.05, antialiased=False)
ax.set_zlim(-1, 1)
ax.zaxis.set_major_locator(LinearLocator(10))
ax.zaxis.set_major_formatter(FormatStrFormatter('%.02f'))
fig.colorbar(surf, shrink=0.5, aspect=5)
ax.view_init(elev=25., azim=350.)
plt.show()
Dado que el Hamiltoniano tiene términos de potencia impar en $y$, las soluciones pueden diverger. Para integrarlas utilizaremos la función odeint de scypy. Tomaremos condiciones iniciales al azar con valores alrededor de $0$ pues hay un punto de equilibrio del sistema en $(0,0,0,0)$. A continuación integramos las soluciones y vemos si estan acotadas.
def soluciones_acotadas(tf=1000,R=100,ntotal=500,xini_max=1.):
trayectorias_acotadas=[]
energias_acotadas=[]
condiciones_inicales=[]
tiempos=np.linspace(0,tf,1000)
for i in xrange(ntotal):
xini=pl.rand(4)*(2*xini_max)-xini_max #tomamos un intervalo simetrico centrado en 0
sol=spi.odeint(ec_mov,xini,tiempos)
if max(norma(sol))<R:
trayectorias_acotadas.append(sol)
energias_acotadas.append(funcion_energia(xini))
condiciones_inicales.append(xini)
return tiempos,condiciones_inicales,trayectorias_acotadas,energias_acotadas
t,x0,sols,energias=soluciones_acotadas()
if len(sols)!=0:
print "Encontramos soluciones acotadas. Las energias correspondienes son:"
print energias
print "\n Los valores extremos son:", min(energias),max(energias)
else:
print "No hay trayectorias acotadas"
Ahora, graficamos las trayectorias encontradas en el espacio de posiciones.
for i in xrange(len(sols)):
plt.plot(sols[i][:,0],sols[i][:,1])
#plt.xlim([0,10])
plt.show()
Observamos que las energías que encontramos son todas positivas y que estas estan confinadas a una region triangular. En la gráfica del potencial, esta parece ser relativamente plana.
Sabemos que para la condición inicial $(0,0,0,0)$ hay un punto de equilibrio, por lo que la trayectoria también esta acotada. Podemos inferir que la mínima energía en la que hay trayectorias acotadas es 0. A partir de la gráfica del potencial, podemos intuir que es posible que existan trayectorias cerradas para energías mayores, pero estas dependen de las condiciones iniciales y deberían encontrarse alrededor del origen, dentro de esta zona.
b) Ahora integraremos las ecuaciones de movimiento con condiciones iniciales tales que $x=0$ y $E<\frac{1}{6}$. Dado $(y,p_y)$ y $E$ se tiene automáticamente un valor para $p_x$. Programaremos una función que haga esto.
def px_inicial(y,py,E):
if E>1/6.:
print "La energía debe ser menor que 1/6."
return np.nan
else:
px=np.sqrt(2*(E-0.5*py**2-0.5*y**2+(1/3.)*y**3))
return px #si el radical es negativo regresa un NaN
Ahora graficamos una trayectoria para ejemplificarlo:
E=0.1
x=0.
y,py=0.3,0.
px=px_inicial(y,py,E)
xini=[x,y,px,py]
tiempos=np.linspace(0,100,1000)
sol=spi.odeint(ec_mov,xini,tiempos)
plt.plot(sol[:,0],sol[:,1])
plt.plot(xini[0],xini[1],marker='o',markersize=10,color="red")
plt.plot(sol[-1,0],sol[-1,1],marker='^',markersize=10,color="red")
plt.grid()
plt.show()
La trayectoria esta acotada entro del intervalo de tiempo en el que se graficó.
c)Ahora, programamos una función que nos de el valor de $p_y$ de manera angáloga al inciso anterior.
def py_inicial(y,px,E):
if E>1/6.:
print "La energía debe ser menor que 1/6."
return "Error"
elif 2*(E-0.5*px**2-0.5*y**2+(1/3.)*y**3)<0:
return "Error"
else:
py=np.sqrt(2*(E-0.5*px**2-0.5*y**2+(1/3.)*y**3))
return py
Ahora programamos una función que entuentre los puntos donde $x=0$ de manera similar a lo que se hizo en la tarea examen anterior.
def trayectoria_con_ceros(inicio,N=100,mi_dt=0.1,error=1e-5):
"""
Realiza integracion de la trayectoria mediante odeint. Si encuentra un cruce con x=0, calcula el punto donde esto
ocurre y continua la integracion.
"""
tiempos=np.array([0.])
trayectoria = np.zeros((N,len(inicio)))
trayectoria[0,:] = np.array(inicio)
raices=np.zeros((0,4))
contador_cruces=0
for i in xrange(N-1):
ti=tiempos[-1]
t_local=np.linspace(ti,ti+mi_dt,N)
sol_dt = spi.odeint(ec_mov,trayectoria[i,:],t_local[1:N])
tiempos=np.concatenate((tiempos,np.array([t_local[N-1]])))
#ahora buscamos si hubo algun cruce
if sol_dt[-1,0]*sol_dt[0,0] <= 0: #si hay cruce con el plano izq
contador_cruces+=1
punto_de_cruce=cruce_con_cero(trayectoria[i,:],mi_dt,N,error)
raices=np.vstack((raices,punto_de_cruce))
#agregamos el punto de cruce a nuestra lista de raices
trayectoria[i+1,:] = sol_dt[-1,:] #tomamos solo el punto final de la integración en la trayectoria
return contador_cruces,raices,tiempos,trayectoria
def cruce_con_cero(punto_actual,mi_dt,N=100,error=1e-5):
dt_busqueda = mi_dt/2
y_izq = punto_actual #punto inicial del proceso de integracion actual
t=np.linspace(0,dt_busqueda,50)
sol = spi.odeint(ec_mov,y_izq,t)
y_med = sol[-1,:] #toma el ultimo punto de la integracion como el valor de la solucion en el punto medio
n=0
while abs(y_izq[0])>error and n<=N:
if y_izq[0]*y_med[0] <= 0:#hubo un cruce
reinicio = y_izq
else:
reinicio = y_med
dt_busqueda = dt_busqueda/2
t=np.linspace(0,dt_busqueda,50)
sol = spi.odeint(ec_mov,reinicio,t)
y_izq = reinicio
y_med = sol[-1,:]
n+=1
return y_izq
Ahora tomaremos una serie de trayectorias con distintos valores iniciales, tales que $x=0$, $y \in[0,0.1]$, $p_x \in[-0.15,0.15]$ y la energía es uno de los siguientes valores: $E=0.01,0.03,0.1,0.12,0.15,0.166$
lista_energias=[0.01,0.03,0.1,0.12,0.15,0.166]
fig=plt.figure(figsize=(15,50))
fig.suptitle('Simulaciones con Condiciones Iniciales Aleatorias',fontsize=18)
#plt.subplots_adjust(wspace=0.2)
plt.subplots_adjust(hspace=0.2)
for i in xrange(15):
E=choice(lista_energias) #elige una energía al azar de la lista
x=0.
y=pl.rand()*0.1
px=pl.rand()*.3-.15
py=py_inicial(y,px,E)
xini=[x,y,px,py]
cont,roots,tiempos,sol=trayectoria_con_ceros(xini,N=500)
ax=fig.add_subplot(8,2,i+1) #haremos un arreglo de 8x2 subplots
ax.set_title('xini=[{},{},{},{}] y E={} \n'.format(x,y,px,py,E),fontsize = 10)
plt.plot(sol[:,0],sol[:,1])
#plt.plot(sol[0,0],sol[0,1],marker='o',color="red")
ax.plot(roots[:,0],roots[:,1],'ro')
plt.ylabel('$y$',fontsize=20)
plt.xlabel('$x$',fontsize=20)
plt.grid()
plt.show()
Dado que utilizamos parámetros aleatorios para constuir la condición inicial, puede ser que la condición de la energía dada a la función py_inicial no pueda cumplirse.
d)Finalmente, tomamos condiciones iniciales de manera similar al inciso anterior, que sean compatibles con la condición de energía pedida y graficamos el espacio fase $(y,p_y)$ para cada una. Sin embargo, con el fin de tener una muestra representativa del espacio fase, tomaremos varias condiciones iniciales para cada valor de la energía. Además, para cada energía y su respectiva condición inicial graficaremos las raices encontradas,i.e, los puntos donde $x=0$.
De este modo, para cada energía se tendra la trayectoria graficada del lado izquierdo y las raices del lado derecho.
lista_energias=[0.01,0.03,0.1,0.12,0.15,0.166]
#plt.subplots_adjust(wspace=0.2)
#plt.subplots_adjust(hspace=0.1)
for i in xrange(len(lista_energias)):
E=lista_energias[i]
x=0.
fig=plt.figure(figsize=(15,7))
fig.suptitle('Energia={}'.format(E),fontsize=18)
ax1=fig.add_subplot(1,2,1) #haremos un arreglo de 2 subplots
ax2=fig.add_subplot(1,2,2)
for j in xrange(10):
py="Error"
while py=="Error":
y=pl.rand()*0.1
px=pl.rand()*0.3-.15
py=py_inicial(y,px,E)
xini=[x,y,px,py]
t=np.linspace(0,100,1000)
sol=spi.odeint(ec_mov,xini,t)
ax1.plot(sol[:,1],sol[:,3])
ax1.set_ylabel('$p_y$',fontsize=20)
ax1.set_xlabel('$y$',fontsize=20)
cont,roots,tiempos,tray=trayectoria_con_ceros(xini,N=4000)
ax2.plot(roots[:,1],roots[:,3],'ro')
ax2.set_ylabel('$p_y$',fontsize=20)
ax2.set_xlabel('$y$',fontsize=20)
plt.show()
Lo primero que observamos en el espacio fase es que las trayectorias comienzan siendo circulares y se van deformando horizontalmente para formar una especie de "trompo" o "huevo". A si mismo, observamos que al aumentar la energía, estas comienzan a invadir el centro de la gráfica. Es decir, llenan una mayor cantidad del espacio fase que previamente era inaccessible a las energías. Por otro lado, observamos que las trayectorias en el centro van perdiendo su simetría circular y adiquieren patrones más irregulares. Aún así, todas las orbitas estan acotadas dentro del intervalo de tiempo utilizado.
En cuanto a las raices podemos observar que estas también tienen acceso a más regiones del espacio conforme aumenta la energía. Cabe mencionar que para las primeras 2 energías se observa que hay raices del sistema donde no pasan trayectorias, esto puede deberse a que realizamos procesos de integración separados para cada uno con tiempos distintos, por lo que las soluciones no son completamente identicas.
Punto Extra¶
Realice un video tipo time-lapse que muestre la evolución de los mapeos para cada uno de los problemas. Puede agregar mapeos intermedios para mejorar la visualización del video.